x = sym('x');
y = sym('y');
z = sym('z');
A = [6 0 0; 0 1 0; 0 0 5.7];
R0 = [-0.2; 0; -0.2];

u0 = sym('u0');
u1 = sym('u1');
u2 = sym('u2');
u3 = sym('u3');

f0 = sym('f0');
f1 = sym('f1');
f2 = sym('f2');
f3 = sym('f3');

u0 = [x; y; z];
u1 = -A * u0 - R0 - R1(u0) - R2(u0, u0);
u2 = -A * u1 - R1(u1) - R2(u1, u0) - R2(u0, u1);
u3 = -A * u2 - R1(u2) - R2(u2, u0) - 2 * R2(u1, u1) - R2(u0, u2);

f0 = R0 + R1(u0) + R2(u0, u0);
f1 = R1(u1) + R2(u1, u0) + R2(u0, u1);
f2 = R1(u2) + R2(u2, u0) + 2 * R2(u1, u1) + R2(u0, u2);
f3 = R1(u3) + R2(u3, u0) + 3 * R2(u2, u1) + 3 * R2(u1, u2) + R2(u0, u3);
